function MVS = selectMeanVelocity(MV,IDV,d1,d2)
%selectMeanVelocity - Mean Velocity of selected reflectors
%
%   MVS = selectMeanVelocity(MV,IDV,d1,d2)
%
% A velocity time series MVS, in terms of absolute velocity, is computed 
% for each time in MV(:,1,1) in the range from d1 to d2 and taking the 
% mean from the targets whose indices are in the index set IDV.
% If d1 (d2) is undefined or empty, d1 = MV(1,1,1) (d2 = MV(end,1,1)) is
% used.
%
% See also extraTarget2cell, cell2posArray, velArray.

if (nargin < 4)||isempty(d2)
    d2 = MV(end,1,1);
end
if (nargin < 3)||isempty(d1)
    d1 = MV(1,1,1);
end

I1 = find(MV(:,1,1) == d1);
I2 = find(MV(:,1,1) == d2);

MVR = MV(I1:I2,:,:);
MVR = MVR(:,:,IDV');
[nr,~,nt] = size(MVR);

MVS = zeros(nr,2);
MVS(:,1) = MVR(:,1,1);

for k = 1:nr
    MVRk = MVR(k,5:7,:);
    absVk = zeros(nt,1); 
    for kk = 1:nt
        absVk(kk) = norm(MVRk(:,:,kk),2);
    end
    Vk = nanmean(absVk);
    MVS(k,2) = Vk;
end

% Search of outliers
MVm = nanmedian(MVS(:,2));
max(MVS(:,2));
Im = find(MVS(:,2) >= 100*MVm);
if ~isempty(Im)
    nm = numel(Im);
    for k = 1:nm
        Imk = Im(k);
        n1k = max(1,Imk-50);
        n2k = min(nr,Imk+50);
        MVSk = MVS(n1k:n2k,2);
        Ink = find(MVSk >= 100*MVm);
        nek = setdiff((1:(n2k-n1k+1))',Ink);
        MVSek = MVSk(nek);
        MVS(Imk,2) = nanmean(MVSek);
    end
end